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SECTION  I 
INTRODUCTION 


For  the  past  several  years,  both  the  Navy  and  industry  have 
expended  much  effort  in  the  area  of  optimum  beamforming.  Impetus  for  this 
work  has  come  from  results  in  underwater  acoustics  and  signal  processing 
areas.  As  various  underwater  acoustics  programs  developed,  it  was  learned 
that  many  of  the  noise  fields  encountered  by  current  acoustic  systems  were 
neither  spatially  random  nor  isotropic,  and  in  many  cases,  assumptions  of 
spatially  random  or  isotropic  noise  fields  did  not  even  remotely  approximate 
the  actual  noise  environment.  As  these  facts  developed,  it  became  clear  that 
optimum  beamforming  offered  potential  gains  that  were  too  significant  to  be 
ignored. 

The  goal  of  optimum  beamforming  is  to  form  a beam  in  the 
desired  direction  while  minimizing  the  interfering  noise  arriving  at  the  array 
from  other  directions.  In  those  cases  where  there  are  spatially  discrete  in- 
terfering noise  sources,  optimum  beamforming  designs  a beam  that  has  nulls 
directed  toward  those  interfering  sources.  The  depth  and  directivity  of  the 
nulls  will  determine  the  amount  of  gain  that  is  actually  achieved  against  a 
noise  field  of  this  kind.  If  the  noise  field  has  no  spatially  discrete  compo- 
nents, the  optimum  beamformer  gives  the  same  performance  as  the  conven- 
tional time-shift-and-sum  beamformer.  An  optimum  beamformer  never 
fails  to  at  least  match  the  performance  of  a conventional  beamformer  since 
the  optimum  beamformer  always  samples  the  noise  field  against  which  it 
must  process  and  uses  this  information  in  the  design  of  its  beam.  Thus,  many 
of  the  arbitrary  qualities  associated  with  time-shift-and-sum  beamforming  are 
not  present  in  optimum  beamforming. 
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Almost  concurrent  with  the  developments  in  underwater  acous- 
tics  were  significant  developments  in  digital  signal  processing  techniques. 

One  of  the  most  significant  developments  was  the  fast  Fourier  transform  (FFT) 
algorithm.  The  FFT  development  prompted  the  evolution  of  a whole  new 
series  of  algorithms  for  the  digital  implementation  of  optimum  beamforming. 

In  the  early  stages  of  Navy  application,  all  optimum  beam- 
forming was  done  in  the  time  domain.  The  time-domain  approach  was  severely 
limited  by  the  large  number  of  arithmetic  operations  required  to  implement  the 
digital  filters.  However,  the  advent  of  the  FFT  opened  new  avenues  for  ap- 
proaching the  problem.  With  the  capability  for  doing  frequency-domain  pro- 
cessing, a series  of  new  algorithms  were  developed.  These  algorithms  were 
much  lqore  reasonable  in  that  they  greatly  reduced  the  number  of  arithmetic 
operationsN^jequired  to  implement  digital  filters. 

This  report  w-iH  present*the  theory  behind  the  frequency-domain 
approach  to  optimum  beamforming  as  implemented  with  Wiener  multichannel 
filters,  along  with  algorithms  for  actually  designing  the  filters.  There  are 
other  approaches  to  optimum  beamforming  that  could  be  taken  with  comparable 
success.  However,  Wiener's  approach  is  the  optimum  under  the  least-mean- 
square  criterion  for  signal  extraction.  For  this  reason,  Texas  Instruments 
has  chosen  to  use  Wiener's  approach  for  initial  evaluation  of  optimum  beam- 
forming  before  considering  other  approaches  that  approximate  optimum.  iy 

Another  subject  proven  valuable  in  studying  the  ability  of  mul- 
tichannel filters  to  perform  is  the  frequency-wavenumber  spectrum.  The  \ 
frequency-wavenumber  spectrum  is  used  to  characterize  the  acoustic  noise 
field  that  the  filters  must  process  against.  Since  this  technique  provides  in- 
formation essential  to  the  understanding  of  optimum  beamforming  performance, 
this  report  will  include  a discussion  of  frequency-wavenumber  spectra. 

Finally,  to  enable  comparison  of  optimum  beamforming  with 
time-shift-and-sum  beamforming,  array  gain  techniques  will  be  discussed. 
Along  with  the  actual  array  gain  calculation,  the  type  of  comparisons  that 
can  be  performed  are  considered. 


SECTION  II 

WIENER  MULTICHANNEL  FREQUENCY- DOMAIN  FILTERING 


To  provide  the  necessary  background  for  understanding  Wiener 
filtering  and  the  problems  encountered  when  actually  computing  the  filters,  a 
complete  theoretical  discussion  of  Wiener  filtering  will  be  undertaken.  A 
derivation  will  be  given  which  provides  the  exact  solution  for  the  Wiener  mul- 
tichannel filters  in  the  frequency  domain.  Since  digital  filtering  on  a computer 
requires  that  filters  be  limited  to  a certain  total  length,  it  will,  in  general, 
be  impossible  to  apply  these  frequency  filters  exactly  because  their  inverse 
Fourier  transform  ordinarily  will  be  of  infinite  length  in  time.  The  derivation 
is  given  because  it  is  needed  to  obtain  an  understanding  of  multichannel  filter- 
ing in  the  frequency  domain. 

The  performance  criterion  in  Wiener  linear  multichannel  fil- 
tering is  to  minimize  the  difference  e(t)  between  the  signal  s(t)  and  the  output 
g(t)  of  a Wiener  filter  set  ja^r)  | j = 1,  2,  ...»  n;  - « < t < ooj  by  minimizing 
the  mean-square -error  I given  by 


I 


lim 
T -•  00 


g(t)  - s(t) 


2 

dt 


where  g(t)  is  the  output  of  the  multichannel  filter  set  obtained  by  convolution 
filtering  of  each  channel  and  subsequent  summation  of  the  filtered  time  traces: 


g(t) 


a^-r)  L(t-t)  dT 
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square -error,  we  obtain 
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(t-x)dtSdrdx 


(t)  between  the  k1*1  channels  at  a 


The  correlation  function 


The  subscript  s will  denote  correlation  with  the  signal  function  s(t) 


The  mean  square  error  I is  then  simplified  to 
n n f°° 

1 23  aj'T>  vt  + *> dT  <pw 

j=i  k*i  J-°=  LJ-°°  J J KJ 


(x)  dx 


r 00 

n r 

•2E  J.„ajM 

j=l 


(T)dT  + cpss(0) 


,» . r , 

J L J 


To  see  how  changing  the  frequency  filter  weights 
-i2rrf  r 

(T)  e dT  affects  the  value  of  the  mean  square  error, 


recall  that  the  Fourier  transform  and  its  inverse  transform  are  related  in 
the  following  manner: 


<«)ei2,,£,d.  ei2n£T  df  = a.(T) 


I ( , 

J -CO  J -CO  J 


Hence,  an  alternative  form  for  the  impulse  response  a.(T)  of  the  frequency 

filter  A .(f)  is  given  by 

J < ® 


i.  (T)  = f A.  (f)  e'?'nfT  df 

J J.»  J 


Recall  also  that  convolution  in  the  time  domain  or  the  frequency 
domain  results  in  multiplication  in  the  frequency  domain  or  the  time  domain 
when  a Fourier  transform  is  performed.  Here,  the  following  result  will  be  used; 


m 00  -00 

jL  [ L ax (T  ■ x)  dT 


-i2nfx,  . . ... 

e dx  = A.(f)  A (f) 

J k 


so  that 


/ /*. 


(T)a  (T+x)dT  e + l2nfxdx  * A .(f)  A,  ( -f)  = A .(f)  A*  (f) 
k J k J k 


where  the  asterisk  denotes  complex  conjugate. 


D 


II-4 


science  services  division 


From  the  preceding  results, 


Since 


n n 


Aj(f)A*  (f)  e'i2TTfxdf 


c^j(x) 


dx 


cp-fr)  dT  + cp  (0) 
JS  ss  ' ' 


-i2rrfx 


V 


dx 


df 


n rc 

-zV  I 


■fsf  j-.V£,l  l_VT)  =i2"£TdT 


df  + tp  (0) 
SS 


n 

* A.(f)  A*(f>  §*k  (f)  - 2 52  A.(f)  is  .,f, 

j = 1 


df 


/•  00  QQ 

jL  AJm  vf>  * - /_.  Aj,-f)  v-£>  * * /./* m *.* 


(f)  df 


it  follows  that 


n n 


j*l  k=l  J Jk 


Aj(f)  »gj(£)  + A*  (f)  **.(f) 


df 


+ *ss<°> 


Thus,  the  quantity  inside  the  braces  is  seen  to  be  a real  function  of  frequency. 
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To  find  the  filters  A^(f)  which  minimize  the  mean-square-error, 
it  is  necessary  to  set  in  the  partial  derivatives  of  the  integrand  with  respect  to 
the  real  and  imaginary  parts  of  A.(f)  = 0 for  each  channel: 


n 


k=l 


Vf>  V£) 


5jt 

A,  (f) 
k 


(f) 

jk 


$ .(f)  - $*.(f) 
sj  sj 


A*  (f)  §*,(f)  - AJf)  *.Jf)  - $_;(f)  + $*,(f) 


k=l 


jk' 


jk' 


sj 


sj 


= 0 


Hence,  the  optimum  filters  A.(f)  must  be  such  that 

J 


£ Vf>  V(£)  • 

k=l 

Therefore,  it  is  necessary  to  solve  the  matrix  equation 


■*!!'«  ?U<4>  • ' ' *!„<*>' 

Aj(f)‘ 

§22<f>  ' ■ ' *2„<f> 

A2(f) 

’25(f) 

Vf>_ 

to  obtain  the  Wiener  optimum  multichannel  frequency  filter  setjA.(f)}. 

J 

This  matrix  equation  forms  a system  of  only  n equations  in 


•1 


0 

0 

D 


the  unknowns  A.(f)  but  must  be  solved  for  all  of  the  frequencies  f in  the  band 
J 

of  interest. 


Assume  that  a different  frequency  filter  set  A.(f)  + A.(f)  is  applied  to  the  par- 

J J 

ticular  set  of  sensor  outputs  for  which  the  correlation  transforms  were  com- 
puted. The  mean-square-error  I will  then  be  given  by 


n 

= L K,f)  - § b( 


bs(f)  - Xj  [A.(f)  ^sj(0  + A*(f)  fyf)' 


E h 

j=l  L J 


n n 


(f)  $ (f)  + A 


* * 1 

j(f) 


+i?  S h<£)  + Ai(f)]  [v£)  + Ak(0]  *jk(1i } 
= L K(£)  ■§  h(f)  v°  + a>  •>' 


n 

-E[- 

j=l  L J 


(f)  § .(f)  + A*(f)  $*.(f) 
J SJ 


n n 

E E k 

j=i  k= i L J 


.(f)  A*(f)  $Jk(f)  + A.(f)  A*(f)  $;k(f) 


+ A*(f)  Ak(f)  $.k(f)  + A.(f)  A*(f)  $*k(f)lldf 
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I 


I 


0 

iD 

C 


m 


L la.,  i i- 


Since 


I 

E vf)  ’jk10  * *j. 

k=  l 

An 

and 

I 

11 

> A (f)  $ (f)  = $.  (f) 

k jk  is' 

«.*• 

k = 1 

n n 

n 

• e 

E E 

vf)  vf>  vf)  - E 

j=  1 k = 1 

j = i 

i * 

and 

n n 

n 

* * 

E E 

Aj(f)  A'(fl  »jk(l,  = Y, 

j=  1 k = 1 

j = l 

and 

n n 

n 

l • 

E E 

j=  1 k = 1 

V°  AkS(£)  ,*k(f)  = £ 

j=l 

A.(f)  i .(f) 
J sj 


the  mean-square -error  I is  given  by 


I = 


/“  n n 

*-lfl  ■ E A"<f>  <> + 32  E va  v« 

- *>  l J=I  j = 1 k = 1 


df 


f *ss(f)  -32  \(f)  v° 

00  j = 1 


df 


r n 


+ lim 


T -e 


-T 


2'T  / tr  / VT,fj 


T2 


6 (t)  f.(t  - T)  dT 


U=  1 


dt 


— 
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I 

! 

T 


i. 


Since  the  bottom  expression  is  positive,  the  mean-square-error  I is  a mini- 


mum if  and  only  if  6^(T)  = 0 for  all  values  of  j and  T • This  condition  is  equiv- 


alent to  the  condition  that  A.(f)  = 0 for  all  values  of  j and  T.  Therefore,  the 

Wiener  optimum  frequency  filter  set  |A.(f)J  is  achieved  by  satisfying  the  set 
of  equations 


n 

Z Vf>  Vf)  - Vf) 

k = 1 


Then  the  mean-square-error  I is  given  by 


n n 


I 


fl 

ii 

o 

H Ic 


i. 


: B 

i i 

Ii 

! i L 


0 
E 

in  b 


£ £ Ml 

j = 1 k = i L J 


-l 


n n r- 

=£  £ M 

j=  1 k= 1 L J 


-1 


n a i 

S £ p>l 

j= 1 k= 1 L J 


*.  (f)  § ,(f) 

js  sk 


Of)  *> 


<s(f)  $sj<f) 


n 


The  quantity 


V»  #„« 

j=l 


V*  A .(f)  $ .(f) 

Z—f  J SJ 

j = l 


is  equal  to  its  complex  conjugate  and  is,  therefore,  real.  It  is  then  reason- 
able to  define  a mean-square-error  density  1(f)  as  follows: 


1(f)  = $ (f) 

s s 


£ £ isj,£,[v,I)l 

j = 1 k= 1 L J 


-1 


*ks(f) 


n n 


= * Q(f) 

s s 


£ £ ♦.jw  Ml 

j = 1 k= 1 L J 


-1 


\JC) 


* (f> 

ss 
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The  quantity 


E E 5sj,f)  k<f)]  * fk9(f) 

j=  1 k= 1 L J 


is  known  as  the  multichannel  squared  coherence  of  the  signal  s(t)  with  the 
sensor  output  set  Jf.(t)  | j = 1,  2, nj  at  the  frequency  f.  It  is  a mea- 

sure of  the  relative  success  of  a Wiener  optimum  multichannel  filter  set 
in  extracting  a signal  s(t)  buried  in  noise.  If  it  is  equal  to  1,  the 
mean-square-error  density  1(f)  at  the  frequency  f is  0.  If  it  is  0,  the  mean- 
square-error  density  1(f)  is  $gg(f)»  the  total  signal  power  density.  In  general, 
it  lies  somewhere  between  these  two  extremes  so  that  the  relative  error 
equal  to  1 minus  the  multichannel  squared  coherence. 


J 


D 

0 

D 

0 

0 

0 


The  analysis  presented  in  Section  II  is  an  analysis  of  contin- 
uous functions.  However,  it  is  not  possible  to  use  continuous  functions  when 
working  with  real  data.  In  this  case,  discrete  values  are  used  to  represent 
the  continuous  functions  and  numerical  methods  are  used  to  approximate 
operations  such  as  integration,  etc.  This  section  will  be  devoted  to  the 
numerical  calculation  of  the  crosspower  matrix  and  the  solution  of  the  matrix 
equation  to  obtain  the  Wiener  filters. 


A.  FAST  FOURIER  TRANSFORM  (FFT) 


D 


0 

0 

0 

I 

0 

0 

1 


Due  to  the  savings  in  computation  time,  all  processing  is  done 
in  the  frequency  domain.  Therefore,  a form  of  the  FFT  is  used  to  transform 
all  time  data  to  the  frequency  domain. 

Suppose  a time  series  having  N = 2n  samples  is  divided  into 
two  series,  each  of  which  has  only  half  as  many  points.  Denote  the  original 
series  y X = (Xq,  X^,  X£,  X3,  ....  Xn)  and  the  other  two  by  Y = (Yq,  Yj, 

♦ • • » j 2)  and  Z = (Zq,  Zj,  ...  respectively.  The  two  smaller  series 

are  related  to  the  original  by 


2k 


k - 0,  1,  2,  ...  , (N/2)  - 1 


Z = X 


2k+  1 


If  the  discrete  Fourier  transform  of  the  time  series  XQ,  Xj, 

...»  X , is  defined  by 
n - 1 7 

Ar  = £ Xk  e-2"i'X/N  r = 0,  1 N - 1 

k = 0 


m-i 
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where  Ar  is  the  rth  coefficient,  Xk  is  the  kth  time  sample  and  j 
the  discrete  transforms  of  Y and  Z are  given  bv 


then 


4rrjrk/N 


4rrjrk/N 


given  in  terms  of  odd  and 


even  numbered  points  as 


4trjrk/N 


For  values  of  r greater  than  N/2,  the  discrete  Fourier  trans 
form  for  and  repeat  periodically  the  values  taken  on  when  r < N/2. 
Therefore,  substituting  r+  (N/2)  for  r in  the  above  equation  for  A gives 
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fin 


OSr<  N/2 


A = 
r 


B + e 
r 


nj.r/ 


Ar+N/2 


e 


2njr/N 


0 s r < N/2 


The  above  process  illustrates  the  reduction  that  can  be  achieved.  By  dividing 
the  sequence  X down  until  the  two-point  transform  is  reached  and  then  using 
the  above  relations,  the  Fourier  transform  of  X can  be  achieved  in  approxi- 
mately 1/2  N • log2  N complex  multiplications  and  N • log2  N additions.  The 
above  process  results  in  a scrambled  Fourier  transform,  making  one  addi- 
tional step  necessary  to  get  the  output  in  the  desired  format. 

B.  CROSSPOWER  MATRIX 

The  crosspower  matrix  is  the  set  of  all  cross-  and  auto-power 
spectral  densities  for  the  various  channels.  Thus,  it  is  necessary  to  compute 
each  of  the  spectral  densities  to  form  the  matrix.  The  method  outlined  by 
Richard  A.  Haubrich  is  a computational  procedure  that  is  particularly  efficient 
for  computing  power  spectra.  * This  method  is  outlined  below. 

Let  Xj(nAt)  be  the  nth  sample  on  channel  j,  where  n = 1,  2,  ....  N. 
Subdivide  X into  p subsets  and  denote  the  kth  subset  by  X.,  (mAt),  where  m = 1, 
2,  ....  (N/p).  Here  it  is  assumed  that  N is  divisible  by  p.  Finally,  denote 
the  Fourier  transform  of  Xjk(mAt)  by  X (f).  The  average  crosspower  spectral 
estimate  between  channels  i and  j is  then  given  by 

P 

Vfl  * ; £ xik«>  x> 

k=  1 J 


Haubrich,  Richard  A.  , 1965:  Earth  Noise,  5 to  500  millicycles  per  second 
J.G.R.  v.  70,  n.  6,  15  Mar.  , p.  H15-1427. 
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where  X^f)  is  the  complex  conjugate  of  X^tf).  Similarly,  the  autopower  for 
the  i^  channel  is  given  by 


*»<« 


Jr 


X.,  (f)  X.,  (f) 
p i—d  lk  jk'  ' 

k = 1 


These  relations  are  equivalent  to  methods  given  by  Blackman  and  Tukey,*  and 

relations  are  equivalent  to  forming  the  correlation  function  Rk.(T)  for  the  kth 
k 

segment  where  R_(T)  is  given  by 


M-|t, 

lij(T)  = M~  X.k(mAt)  Xjk(mAt  + )t|) 

m = 1 


and  then  taking  the  Fourier  transform  of  the  product  h(-r)  RK.(t),  where  1i(t) 
is  a time  window.  For  this  process,  Mt)  is  given  by 


Mt>  - i - i . 


= 0, 


T=0,  ±1,  ...,  ±M 


t = ± (M  + 1),  ±(M  + 2), 


k 

Note  that  for  each  R.  .(t)  formed,  a spectral  estimate  §.  .(f)  is  obtained.  The 

^ J i j 

spectral  estimates  are  then  averaged,  thereby  increasing  the  reliability  and 
stability  of  the  spectral  estimate.  The  final  estimate  .(f)  is  equivalent  for 
both  processes.  The  full  value  of  the  FFT  with  Haubrich's  method  is  easily 
seen  in  terms  of  increased  efficiency  over  the  Blackman  and  Tukey  method. 


Blackman,  R.  B.  , and  J.  W.  Tukey,  1958:  The  Measurement  of  Power 
Spectra,  Dover  Publications,  Inc.,  New  York. 
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C.  SOLUTION  OF  THE  MATRIX  EQUATION 

In  Section  II  it  was  shown  that  the  Wiener  filters  are  obtained 
by  solving  a matrix  equation  of  the  form  AX  = F,  where  A is  a square  Her- 
mitian  matrix  and  X and  F are  both  column  vectors.  The  components  of  the 
column  vector  X are  the  filter  weights;  the  column  vector  F is  the  signal 
model,  and  the  matrix  A is  the  crosspower  matrix  computed  using  the  tech- 
niques previously  discussed.  All  that  remains  is  to  solve  for  X in  terms  of 
the  known  quantities  A and  F.  The  method  currently  used  and  computationally 
more  efficient  than  any  other  technique  used  to  date  is  the  Square-Root  Method. 
This  method  is  approximately  four  times  faster  than  the  Gauss'  Method,  which 
was  the  last  technique  used. 

Since  A is  a square,  Hermitian  matrix,  it  can  be  expressed 
as  the  product  of  two  triangular  matrices  of  which  one  is  the  transpose  of 
the  other.  Thus,  let  A be  given  by 


where 


S11  S12 


S = 0 


S22  S2n 


0 . . . S 


Forming  the  product  S'S  gives  the  elements  of  A in  terms  of  the  elements  of  S. 
The  elements  of  A are  given  by 


= Sn  Su  + S_.  S_.  + ...  + S..  S.„ 

ij  li  lj  2i  2j  n ij* 


a..  = S?.  + S*  + . . . + S2.  , 
ii  li  +i  ii 
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Solving  the  preceding  equations  for  the  S.  . gives 


By  expressing  A as  S S,  solving  the  equation  AX  — F now  becomes  equivalent 
to  solving  the  two  equations 


The  elements  of  K are  determined  by  recurrence  relations  similar  to  those 
given  above.  They  are 


0 

0 


0 


II 


I 


[ 

[ 

[ 

I 

I 

I 

I 


calculations  required  to  find  the  X.. 

l 

This  section  outlined  the  computational  procedures  for  imple- 
menting the  design  of  Wiener  filters  on  a general  purpose  computer.  Many 
factors  influence  how  closely  the  filters  designed  by  the  above  techniques  ap- 
proximate the  theoretical  filters.  However,  by  carefully  applying  the  digital 
techniques,  the  approximation  can  be  made  very  close. 





SECTION  IV 

FREQUENCY  WAVENUMBER  SPECTRA 


A particularly  useful  tool  in  the  study  of  Wiener  filter  perfor- 
mance is  the  frequency-wavenumber  spectrum.  The  frequency-wavenumber 
spectrum  provides  qualitative  information  about  the  noise  field.  This  technique 
provides  a method  for  determining  the  relative  distribution  of  energy  as  a 
function  of  both  frequency  and  wavenumber,  and  the  wavenumber  parameter  is 
directly  relatable  to  spatial  parameters.  Thus,  the  energy  distribution  of 
the  noise  field  can  be  represented  in  terms  of  space  and  frequency.  Section 
IV  presents  a theoretical  discussion  of  the  frequency-wavenumber  spectrum 
and  a description  of  how  the  spectrum  can  be  computed;  other  computational 
techniques  will  not  be  discussed. 

A.  FREQUENCY-WAVENUMBER  POWER-DENSITY  SPECTRUM 

Consider  a plane  wave  of  wavelength  X propagating  in  the  di- 
rection specified  by  the  direction  cosines  (y,  . . . , y ) in  an  n-dimensional 

1 n 

space.  Let  its  time  waveform  at  the  origin  by  given  by 


g(t)  = cos  2n(fQt  + c) 

and  let  kQ  be  the  vector  wavenumber  1/A.  (y^  . . , y^),  which  points  in  the  di- 
rection of  propagation  of  the  plane  wave  and  has  the  magnitude  1/X. 

Since  the  magnitude  of  the  velocity  of  propagation  V is  equal 
to  the  frequency  times  the  wavelength  X and  since  the  perpendicular  distance 
to  any  point  x from  a plane  through  the  origin  with  a normal  specified  by  the 
direction  cosines  (Y1  . . . , Yn>  is  given  by  the  formula 


(Yx  . . . , Yq) 
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the  plane  wave  will  arrive  at  the  position  of  the  point  x after  a time  lag 


•y  ) . x 
n — 


k . x 

= -o  - 


Therefore,  the  equation  of  this  plane  wave  at  position  x is  given 


g(x,  t)  = cos  2rr(f^  t + c - k^.  x) 


If  is  the  vector  displacement  between  the  two  positions  x + y_ 
and  x,  the  space-time  correlation  function  cp(^,  t)  is  the  average  of  the  pro- 
duct of  the  two  functions  g (x,  t)  and  g (x  + t + t)  over  all  values  x of  posi- 
tion and  all  values  T of  time: 


cp  (jr,  t)  = lim  1 

T-*»  _n+l 
. 2 TL 

L-4oo 


T L L 


_ r r r g(*»t ) g (*  + 1 + t 

J-tJ -L  J-L 


) dx  dt 


where  dx  is  vector  shorthand  for  dx, dx  . 

1 n 

Since 


g (x,  t)  = cos  2tt  (fQ  t + c - . x) 

= Re  [ei2"(V  + C ' ^0  ’ £>] 


g (x  + jr,  t + t)  = cos  2rr  (fQ  t + fQ  t + c - k^-x-  k^.^) 


I"  >2TT(f  t + f„T  + c-  k.x-k_.y)l 
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i 


n 


1. 


li 


the  space-time  correlation  function  Vfjr,  t)  reduces  to 


<P  (j£.  t)  = | Re  |^el2TT(f0  T " i^o  * 

after  (n  + l)-fold  integration  over  the  variables  x. x and  T 

1 n 


Let  (t)  and  (t)  denote  the  time  waveforms  of  the  plane 
at  the  positions  x.  and  x^,  respectively.  Then  the  time  correlation  function 


wave 


cp  (t)  = lim 
11  T*. 


2’t  L 


gj  (t)  gt  ( t + T)  dt 


between  the  functions  and  g^  at  a time  lag  r is  also  equal  to 


i Re  ['  ‘2T’(f°T  ' J - i CO.  <£0t  - ^ . z) 


where  jr  is  the  vector  displacement  x^  - x When,  as  in  this  instance,  the 
time  correlation  function  Cfy(T)  between  the  time  waveforms  at  two  positions 
-j  and-^  does  not  depend  upon  their  absolute  locations  - but  only  upon  their 
relative  locations,  together  with  the  delay  time  t — the  wavefield  is  said  to 
be  space -stationary  (in  the  strong  sense).  Any  superposition  of  plane  waves 
is  space-stationary.  For  space -stationary  wavefields,  the  space-time  cor- 
relation function  cp  (£,  t)  j.s  equal  to  the  time  correlation  function  cp.^(T)  be- 
tween any  two  points  x.  and  x ^such  that  £ = x ^ - x . This  fundamental  point 
will  be  a key  building  block  in  a later  discussion. 

0 

0 

I 

K 


The  space-time  to  frequency- wavenumber  Fourier  transform 
5(f,  k)  of  the  space-time  crosscorrelation  function  cp(jr,  T),  called  the  frequency 
wavenumber  power  density  spectrum,  is  defined  to  be 


T L L 


* t”.  f f—f  e 

L-»®  tj'  sl  ■'-L 


-i2n(fT  - k . jr ) 


d^dT 


where  d£  is  vector  shorthand  for  dyj dy^.  ^X>  t)  = £ cos  2tt(]Sq*  X ~ 


f(f.k) 


s(f'fo'  - -v 
2 


+ 6«  + V is  * v 
2 


where  6 is  the  Dirac  delta  function. 

Figure  1 shows  the  frequency-wavenumber  power  density 
spectrum  of  a plane  wave  propagating  in  a 1 -dimensional  space. 


Figure  1.  Frequency- Wavenumber  Power  Density  Spectrum  of  Plane  Wave 
Given  by  Equation  g (x,  t)  = cos  2TT(f0  t + c - kQ  x) 
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jJhawMW 


This  result  is  an  extension  to  n + 1 dimensions  of  the  frequency 
power  density  spectrum 


$(f)  = 


8«-f0>  + 


for  the  sinusoidal  time  function  g(t)  = cos  2n(f^t  + c). 

B.  SPECTRAL,  WINDOW 

When  frequency-wavenumber  spectra  are  computed  using  the 
outputs  g.(t)  of  a set  of  array  elements  specified  by  the  position  vectors 
|xj  | j = i.  2,  ....  m|,  it  is  only  possible  to  obtain  an  estimate  of  the  true 
frequency-wavenumber  spectrum  #(f,  k).  Instead  of  computing  the  integral 

T L L 

H!'~ = iE  / /’/  ^ T>  e"i2T,<fT  - • x,^r 

-T  -L  -L 


it  is  onl<  possible  to  compute  a summation  over  the  vector  displacements  _y 
corresponding  to  the  various  sensor  pairs  of  a spatial  array.  The  set  of 
unique  vector  displacements  between  all  possible  combinations  of  sensor  pairs 
is  called  the  set  of  points  in  correlation  space  corresponding  to  the  array.  If 
several  sensor  pairs  lie  at  the  same  relative  vector  displacements,  they  are 
said  to  be  redundant. 

If  a wavefield  is  composed  of  a superposition  of  plane  waves, 
it  is  space-stationary  in  the  strong  sense;  therefore,  the  space-time  corre- 
lation function  cp(x>  T)  may  be  obtained  by  computing  the  time  correlation 
function 


jt  t-« 


fT 

2T  J gj{t)8^< 


i(t+T)dt 


science  services  division 


between  the  sensor  outputs  g (t)  and  g (t)  of  two  sensors  located  at  x and  x 

J -j  -V 

respectively,  provided  that  the  vector  displacement  x,  - x.  is  equal  to  y. 

V * J 

This  fact  permits  a spatial  sampling  of  the  space-time  correlation  function 
<¥>(x.  T)  over  the  set  of  points  in  correlation  space  corresponding  to  an  array. 

If  the  sensor  outputs  of  an  array  are  not  approximately  consis- 
tent with  a space-stationary  wavefield,  the  frequency-wavenumber  spectrum 
is  not  an  appropriate  analysis  technique.  The  sensor  outputs  are  then  best 
characterized  in  terms  of  the  auto-  and  crosspower  spectra  $.  (f)  between 
sensors. 

A common  method  of  computing  frequency- wavenumber  spectra 
is  to  compute  the  power  spectra 

V(f)  -f_.  ,W2r5rT>«‘I*,fTdT 

and  then  the  summation 

II  »jt -*t> 

j <■ 

to  complete  the  space- time  to  frequency- wavenumber  Fourier  transform.  (If 
several  correlation  space  vectors  x.  - x^  are  redundant,  their  averaged  power 
spectrum  is  used  only  once  in  th’s  summation.) 

This  summation  is  equivalent  to  multiplying  the  Fourier  trans- 
form integrand  by  a sum  of  spatial  Dirac  delta  functions 


E 5(i  - Xm' 
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■ 


where  the  y are  the  various  vector  displacements  comprising  the  correlation 
space  corresponding  to  the  array  used: 


t”  / T / L-  -/L^  t»  e'i2n(,T  ' - ' 11  [ E 6<i  - £„»]  <*<*' 

L-»  -T  -L  -L 


*5  ♦,  fe,-*,, 


J t 


^1*.  • (it  * ij’ 


J l 


Here 


where 


*2  (f)  = l CP(lm>  T)  e"l2TT£T<iT  * MO 

m * 


1*1  X 

T)  = T_«  Jr  J 0 g(x^,  t + t) 


because  we  let 

" 5,  ‘ ^ 

A well-known  result  from  the  theory  of  operational  calculus 
known  as  the  convolution  theorem  states  that  the  complex  Fourier  transform 
of  the  product  of  two  functions  of  a real  variable,  one  of  which  may  be  a 
generalized  function  such  as  the  Dirac  delta  function,  is  the  convolution  of 
the  individual  Fourier  transforms  of  each  of  the  functions: 
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fli 

J £(x)  g(x)  e"l2nyxdx  = Ff(y)  * F 
= Ff(z)  F^(y  - z)  dz 


= y*  Ff(y  - z)  Fg  (z)  dz 


where  F (y)  is  the  Fourier  transform 


f(x)  e-i2"yxdx 


and  F (y)  is  the  Fourier  transform 
© 


/ g(x)  .'“"'"‘dx 

— 00 

The  n-fold  application  of  this  result  implies  that  multiplication 
by  the  quantity  £ 6(x  ‘ corresponds  to  convolution  of  the  true  frequency- 

wavenumber  power  spectrum  #<f,k)  with  the  Fourier  transform  of  £ ^ > 

which  is  equal  to  £ e*2"-  * *m.  The  quantity  W(k)  s £ei2TTii • is  knovS 
as  the  spectral  window  of  the  array  used  in  computing  the  spectrum  $ (f,  k). 

Thus,  the  frequency-wavenumber  spectrum  obtained  will  be 

• a w 

I f *(f’  - ' w -,d-  •/  / #(f.  « ) w(k  - o d^ 
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instead  of  the  true  frequency-wavenumber  spectrum  $(f,k).  If,  at  a partic- 
ular frequency,  all  of  the  true  power  is  concentrated  at  the  wavenumber  k 

—0 

the  frequency- wavenumber  power  density  spectrum  obtained  will  be  a constant 
times  the  wavenumber  spectral  window  W (k)  shifted  by  kQ  in  wavenumber  space. 

Since  the  spectral  window  is  the  sum  of  a finite  number  of  com- 
12TTk  * V 

plex  sinusoids  e - it  cannot  be  the  spike  6(k)  required  to  eliminate  dis- 
tortion of  the  true  frequency-wavenumber  spectrum.  In  fact,  to  achieve  the 
ideal  spectral  window,  i.  e.  , the  Dirac  delta  function  6 (k),  the  correlation 
space  sampling  function  must  be  a function  equal  to  1 at  all  points  of  the  n- 
dimensional  correlation  space  (y  x,  ....  yn).  The  spectral  window  of  any  array 
of  discrete -point  sensors  will  have  rounded  peaks  over  all  of  wavenumber  space. 
The  peak  at  k = £ is  called  the  main  lobe.  All  other  peaks  are  called  sidelobes. 
To  minimize  the  distortion  of  frequency- wavenumber  spectra,  the  spectral 
window  sidelobes  should  be  as  low  as  possible. 

The  same  method  can  be  applied  to  determine  the  angle  of  in- 
cidence for  waves  propagating  across  a 2-dimensional  array.  For  a particular 
frequen'  y,  it  will  be  possible  to  draw  rings  of  constant  apparent  velocity  or 
(equivalently)  of  constant  angle  of  incidence.  The  angle  of  incidence  0 will  be 
determined  similarly  to  the  1- dimensional  case  by  the  formula 

0 = sin*1^  |k|  IVpJ/fJ  = sin'^IVj/IVj) 

Each  ring  will  be  centered  at  the  origin  of  the  wavenumber  spectrum  for  a 
particular  frequency  and  will  be  the  circle  determined  by  taking  a horizontal 
slice  at  the  frequency  f of  the  corresponding  constant- velocity  cone  in  the 
3- dimensional  representation  of  the  frequency- wavenumber  spectrum. 


I 
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c.  COMPUTATION  OF  FREQUENCY- WAVENUMBER  SPECTRA 


The  computation  of  a frequency-wavenumber  power-density 
spectrum  • (f,k)  is  similar  to  the  computation  of  a frequency  power-density 
spectrum  $ (f).  In  the  case  of  a frequency  power-density  spectrum,  a Fourier 
transform  from  time  t to  frequency  f is  made  while,  in  the  case  of  a frequency 
wavenumber  spectrum,  a Fourier  transform  from  space  and  time  (jr,  t)  to 
frequency  and  wavenumber  (f,  k)  is  performed.  In  the  conventional  power- 
density  spectrum  computation,  the  autocorrelation  cp(T)  is  transformed  to 
yield  the  power  density  spectrum 


$ (f) 


-i2rr  fT 

cp  (t)  e dr 


On  the  other  hand,  in  the  frequency-wavenumber  power  density  spectrum 
computation,  the  space-time  correlation  function  qp (_jr,  t)  is  transformed  to 
yield  the  frequency- wavenumber  spectrum 


00  CD  00 


*(f. 


h)  - j J'~'J  V T) 


-i2n  (fT-k  . x.) 


djrdT 


Here,  as  mentioned  earlier,  n indicates  an  n-fold  integration  corresponding 
to  the  number  of  dimensions  of  the  array,  and  d£  is  vector  shorthand  for 
dyi  -~dV 

E&rliQrt  it  was  shown  that  this  definition  of  a frequency-wave- 
number  power  density  spectrum  does  indeed  yield  a physically  meaningful 
representation  of  the  space-time  wavefield.  There,  a frequency-wavenumber 
power  - density  spectrum  f (f,  k)  was  shown  to  transform  a unit-amplitude  plane 
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wave  of  frequency  f and  wavenumber  k into  two  half- strength  impulses  (i.  e. , 
Dirac  delta  functions)  at  (f,k)  and  (-f,  -k)  just  as  a frequency  power -density 
spectrum  transforms  a unit- amplitude  sinusoidal  waveform  of  frequency  f in- 
to two  half- strength  impulses  at  f and  -f. 

The  algorithm  for  the  transformation  of  a space-time  correlation 
function  cp  (y , t)  to  the  corresponding  frequency-wavenumber  power  density  spec- 
trum $(f,k)  consists  of  two  steps. 

The  first  step  is  the  computation  of  the  crosspower  statistics 
of  the  wave  energy  in  the  form  of  interchannel  crosspower  functions  between 
the  outputs  of  the  elements  in  the  sensing  array.  The  techniques  for  computing 
the  crosspower  spectra  is  the  same  as  that  discussed  in  Sec.  III. 

The  second  step  is  the  summation 


£ V ,f)  e 


-i2rr k . (x.  - x,  ) 
- -J  -l 


which  completes  the  transformation. 


Provided  that  energy  travels  unattenuated  across  the  array  in 
the  form  of  plane  waves  with  negligible  refraction,  reflection  or  diffraction, 
the  crosspower  spectra  (f)  corresponding  to  a particular  vector  displace- 
ment will  be  identical.  This  assumption  of  strong  space-stationarity  will  be 
approximately  correct  if  energ>  sources  are  distant  from  the  array  and  the 
medium  in  the  vicinity  of  the  array  is  near-uniform.  The  averaging  of  all 
transforms  ^ (f)  corresponding  to  a particular  vector  displacement  will  tend 
to  reduce  the  effects  of  non- space  stationarity. 
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SECTION  V 
ARRAY  GAIN 


In  the  investigation  of  Wiener  multichannel  filtering,  it  is  ne- 
cessary to  have  some  oasis  for  comparison  between  multichannel  filter  per- 
formance and  the  performance  of  a well-understood  standard.  The  basis  chosen 
is  the  amount  of  array  gain  achieved  compared  to  the  gain  of  a time- shift  - 
and-sum  beamformer.  Time-shift-and-sum  beamforming  is  the  type  of 
beamforming  currently  used  in  a great  many  acoustic  systems  and  is  well 
understood  by  a great  many  people.  Therefore,  comparing  the  amounts  of 
array  gain  obtained  by  Wiener  filtering  to  that  for  time-shift-and-sum  beam- 
forming is  both  realistic  and  meaningful. 

The  array  gain  of  a particular  processing  system  is  defined 
to  be  the  output  signal-to-noise  ratio  divided  by  the  input  signal-to-noise  ratio 
of  a single  sensor.  For  multichannel  filtering  or  time-shift-and-sum  beam- 
forming,  the  output  signal-to-noise  ratio  for  a given  frequency  is  given  by 
the  following  matrix  equation: 


(S/N) 


AH  fSA 

ah  *na 


s 

where  $ is  the  signal  crosspower  spectrum  matrix,  $N  is  the  noise  cross- 
power spectrum  matrix,  A is  a column  matrix  whose  components  are  either 
the  optimum  filter  set  or  the  time-shift-and-sum  filter  set,  and  H denotes 
the  transpose  conjugate.  If  the  matrix  A consists  of  the  time- shift- and  sum 
filter  set,  then  A is  given  by 

iZrrf  T i 
e 

i2rrf  T 0 


i2TTf  T 


where  each  i\iti  = 1,  2,  ...»  n,  is  the  time  delay  required  to  align  an  incom- 
ing plane  wave  in  phase  so  that  coherent  addition  will  occur  when  the  outputs 
of  all  the  sensors  are  summed. 

There  is  a problem  when  it  comes  to  defining  the  input  signal- 
to-noise  ratio.  If  one  sensor  in  the  array  is  chosen  from  which  the  input 
signal-to-noise  ratio  is  determined,  then  the  assumption  has  been  made  that 
the  noise  level  is  the  same  at  all  other  sensors.  This  assumption  is  hard  to 
validate.  For  this  reason,  the  input  signal-to-noise  ratio  is  taken  to  be 


(S/N).n  = 


- t *> 

n j«I  > 


where  $ (f)  is  the  signal  autopower  spectrum  and  $ . (f)  is  the  noise  autopower 

til  ^ 

spectrum  as  measured  by  the  j sensor. 

Using  the  above  definitions  for  input  and  output  signal-to-noise 
ratios  gives  the  array  gain  measurement  G(f) 


G(f)  = 


(S/N) 

t 

(S/N) 


= AH(f)  $S(x)  A(f)  / $S(f) 

AH(f)  §N(f)  A(f)  j JL  £ $N(f) 

j=l  J 

Let  A be  the  matrix  of  optimum  filter  weights  and  let  B be  the 
matrix  of  time- shift-and- sum  filter  weights.  To  compare  the  performance 
of  multichannel  filtering  with  that  of  time-shift-and  sum  beamforming,  the 
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ratio  of  the  array  gain  achieved  by  optimum  beamforming  to  that  achieved  by 
time- shift-and- sum  beamforming  is  formed.  This  ratio,  denoted  by  C(f),  is 
given  by 


C(f)  = 


AH  ftSA 

A § A 

H XN 
A $ A 


H _S_ 
B $ B 

H N 
B $ B 


The  function  C(f)  can  be  used  to  compare  the  performance  of  the 
two  systems  when  trying  to  process  against  the  same  noise  field.  The  main 
advantage  achieved  by  taking  the  ratio  of  the  array  gains  is  that  C(f)  is  not 
dependent  on  the  input  signal-to-noise  ratio.  Thus,  the  outputs  of  the  proces- 
sor can  be  directly  related. 
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